#!/usr/bin/perl -w
use strict;

use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin", "$FindBin::Bin");
use db_parameters;
use ensembl_parameters;
use chromosomes;
use min_max;
use die_error;
use run_cmd;
use DBI;
use align_gaps;
use Data::Dump qw(dump);

print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    table_exon_sequences
	

	Description:
		Obtain the exonic DNA associated with each transcript
	Populates DB tables:
		prot_exon_sequences
		exon_sequences

888888888888888888888888888888888888888888888888888888888888888888888
HEADLINE

use constant GENEID =>0;
use constant CHRMSM =>1;
use constant START_ =>2;
use constant END___ =>3;
use constant STRAND =>4;
use constant CENTRE =>5;
use constant ORDER_ =>5;


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Determining exon structures of coding sequence


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#_________________________________________________________________________________________

#	open_chrm_file

#		retrieve genomic sequence

#_________________________________________________________________________________________
sub open_chrm_file($$)
{
	my ($taxon_data, $chrm_name) = @_;
	#
	#	return null if file not found
	#
	unless (open(DNA, "$taxon_data->[SEQ_DIR]/dna/chromosome.$chrm_name.fa"))
	{	warn
		"Error:\n\tCould not open '$taxon_data->[SEQ_DIR]/dna/chromosome.$chrm_name.fa'\n";
		return undef;
	}
	my $ignore = <DNA>;
	return (*DNA, length($ignore));
}




#_________________________________________________________________________________________

#	retrieve_exon_coordinates

#_________________________________________________________________________________________
sub retrieve_exon_coordinates($$)
{
	my ($dbh, $prot_id) = @_;
	return $dbh->selectall_arrayref(<<"PL/SQLCMD");
	SELECT
		prot_id,
		prot_exons.start,
		prot_exons.finish,
		strand,
        exon_id
	FROM
		taxon.prot_exons NATURAL JOIN
        ens_id JOIN 
		taxon.gene_loci USING (gene_id)
	WHERE
        prot_id = '$prot_id'
	ORDER BY 
		prot_exons.start;
PL/SQLCMD
}	

#_________________________________________________________________________________________

#	retrieve_exon_coordinates

#_________________________________________________________________________________________
sub retrieve_taxon_chromosome($$)
{
	my ($dbh, $prot_id) = @_;
	my $data = $dbh->selectall_arrayref(<<"PL/SQLCMD");
	SELECT
		chromosome,
		species
	FROM
		ens_id NATURAL JOIN
        taxon.gene_loci
    WHERE
        prot_id = '$prot_id'
PL/SQLCMD
	return @{$data->[0]};
}	


#_________________________________________________________________________________________

#	save_prot_coding_structure

#		retrieve exon coding structures for all proteins in a species

#_________________________________________________________________________________________

sub get_exon_sequences($$\%)
{
	my ($dbh, $prot_id, $taxa_data) = @_;
	print STDERR "\t$prot_id Exons\n";
	my ($chrm, $taxon) = retrieve_taxon_chromosome($dbh, $prot_id);
	print STDERR "$taxon: chromosome $chrm:\n";
	my $taxon_data = $taxa_data->{$taxon};

	{
		print STDERR ("\t\tChromosome $chrm: Reading ");
		#
		#	exon coordinates ordered by protein
		#
		my $exon_coordinates = retrieve_exon_coordinates($dbh, $prot_id);
		die "exons for $prot_id not found" unless scalar @$exon_coordinates;

		#
		#	chromosomal DNA sequence
		#
		my ($chromosome_FH,  $file_offset) = open_chrm_file($taxon_data, $chrm);
		die unless $chromosome_FH;

		
		# for last exon of last protein
		push(@$exon_coordinates, ["LASTLAST"]) if @$exon_coordinates;


		my $cnt_seq = 0;
		my $curr_prot_id = '';

		# exon sequence making up current protein
		my @this_prot_exons_seq;
		my $prot_exons_seq;

		
		$dbh->do("COPY taxon.exon_sequences".
					"(exon_id, dna_sequence) FROM STDIN");

		# for each exon
		for my $line(@$exon_coordinates)
		{
			my ($new_prot_id,  $start, $finish, $strand, $exon_id) = @$line;

			if ($new_prot_id eq 'LASTLAST')
			{
				$prot_exons_seq = join ("",  @this_prot_exons_seq);
				last;
			}

			# get DNA sequence of exon
			#my $exon_seq = substr($chromosome_seq,  $start - 1,  $finish - $start + 1);
			my $exon_seq;
			seek($chromosome_FH, $start - 1 + $file_offset, 0);
			read($chromosome_FH, $exon_seq,  $finish - $start + 1);
			$exon_seq = uc $exon_seq;


			#
			#	3' goes the other way: reverse and complement
			#
			if (!$strand)
			{
				print STDERR "Reverse strand";
				print ">$exon_id unreversed\n$exon_seq\n";
				$exon_seq =~ tr/ACGTacgt/TGCAtgca/ unless ($strand);
				$exon_seq = reverse $exon_seq;
				unshift(@this_prot_exons_seq, $exon_seq);
			}
			#
			#	Save 5'
			#
			else
			{
				push(@this_prot_exons_seq, $exon_seq);
			}

			print ">$exon_id\n$exon_seq\n";
		}
		print ">$prot_id\n$prot_exons_seq\n";
		close $chromosome_FH;
	}
}







#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

# Main logic

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


my $prot_id = $ARGV[0];
chomp $prot_id;
#'mono1|371813|371813|SG';
my $dbh = connect_to_panda();

#
#	Get DB parameters
#
my %taxon_data = get_taxon_data;
get_exon_sequences($dbh, $prot_id, %taxon_data);



$dbh->disconnect();


